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CHAPTER  I 
INTRODUCTION 

1.1  Hot  Dry  Rock  Concept 

Increased  demand  for  versatile  energy  supplies  after  the  energy 
crisis  has  led  to  a  corresponding  interest  in  geothermal  energy.  The 
source  is  located  within  the  earth's  crust  and  methods  used  for  the 
production  can  be  considered  geothermal  extraction. 

In  conventional  geothermal  systems  the  convecting  medium  for 
extracting  the  available  thermal  energy  to  surface  conversion 
facilities  is  the  in-situ  geothermal  fluid  (Palen  and  Narasimhan, 
1981) .  Thus  conventional  systems  depend  on  the  location  of  geothermal 
fluid  reservoirs.  Unfortuna'tely ,  they  are  only  represent  a  small 
fraction  of  geothermal  energy  in  this  world. 

The  Hot  Dry  Rock  (HDR)  concept  does  not  require  the  presence  of 
an  in-situ  fluid  reservoir  but  is  dependent  only  on  the  presence  of  a 
high  geothermal  gradient  and  an  adequate  supply  of  a  working  fluid, 
such  as  water.  It  is  based  on  closed-loop  circulation  of  pressurized 
water  through  a  man-made  fracture  system,  created  by  hydraulically 
fracturing  hot  rock  between  two  wellbores.  The  useful  heat  from 
superheated  water  is  recovered  at  the  surface  through  heat  exchangers , 
and  the  cooled  water  is  reinjected  to  recirculate  through  the 
underground  loop  (Los  Alamos  annual  report,  1984). 

1.2  Recent  Development  In  HDR  Program 

The  Hot  Dry  Rock  concept  originated  at  Los  Alamos  National 
Laboratory   in  1970.   The  major  project  effort  during  the  past  several 


years  has  been  to  develop  a  commercial  size  underground  reservoir  by 
hydraulic  fracturing.  The  program  is  now  largely  centered  on  the 
Fenton  Hill  Project  near  Los  Alamos.  The  site  is  located  on  an 
extinct  volcano  in  the  Jemez  Mountains  of  northern  New  Mexico,  USA. 
The  Phase  II  engineering  system  in  Fenton  Hill  has  a  pair  of  wells. 
EE-2  is  the  injection  well  and  EE-3  is  the  production  well  of  the 
system  (Los  Alamos  annual  report,  1984).  Figure  1.1  shows  a  schematic 
drawing  of  the  reservoir. 

3 
From  December  6   to   9,   1983,   a  total  of  21,300  m  water  was 

injected  into  EE-2   in  61  hours.   This  experiment  was  terminated  by 

fatigue   failure  of  hardware  equipment.    The   resulting  rapid  vent 

returned  to  the  surface  about  54%  of  the  water  that  had  been  injected 

and  delivered  thermal  energy  to  the  surface  at  rates  estimated  from 

100  MW   initially  to  an  average  of  30  MW  over  the  3.3  day  period  of 

rapid  venting.   The  volume  and  rate  of  fluid  return  indicated  that  the 

fracture   system  was   tightly  contained  and  well  connected;  the  high 

rate  of  energy  production  indicated  the  heat  was  extracted  effectively 

from  the  fractured  reservoir  rock  back  through  well  EE-2  (Franke  and 

Nunz,    1985).    Despite  the  disappointment  of  hardware  failures 

preventing  satisfactory  connection  between  two  wellbores,  the  result 

at   least   indicated  that  a  thermal  reservoir  extensive  enough  to  be 

commercially  useful  had  been  opened. 

During  May  and  June  1986,  a  one  month  flow  test  of  this  Phase  II 

3 
heat  extraction  loop  was  conducted.   A  total  of  37854  m  (10  million 

gallons)  of  water  at  20  C  was  injected  into  an  243.8  m  (800  feet) 
long  section  at  depths   around  3657.6  m   (12,000  feet),  where  the 


initial   rock  temperature  was  about  240  C.   Under  a  pumping  pressure 

7     2  3 

2.76x10   N/m   (4000  psi) ,  the  injection  rate  was  0.0179  m  /sec  (285 

gpm) .    The  rate  of  fluid  flow  from  the  production  well  increased  with 

3 
time   to  0.0148  m  /sec   (235   gpm).   The  temperature  of  the  produced 

water  rose  to  190  C.   The  rate  of  heat  production  could  be  converted 

to  about  1500  to  2000  KW  of  electrical  power.   The  rate  of  water  loss 

was   initially  very  high,  in  part  because  of  leakage  through  damaged 

casing   in  the  production  well,  but  it  decreased  with  time  to  a  final 

value   of  26%.   Much  of  this  "loss"  was  water  stored  temporarily  in 

fractures  that  were  outside  of  the  circulation  paths  and  a  large 

fraction  of   it  was   recovered  later  when  the  system  was  vented  (Los 

Alamos  annual  report,  1985). 

In  December,   1987,   experiment  No.  2074  lasted  eight  days.   The 

7    2 
pumping  pressure  was  2.2x10  N/m   (3200  psi)  and  the  outlet  pressure 

was  1.7x10  N/m  (250  psi).  The  input  flowrate  was  96  gpm  and  the 
output  flowrate  rose  from  0  gpm  at  the  start  to  60  gpm  at  the  end  of 
the  test.  The  input  flow  loss  is  believed  to  be  due  to  far  field 
leakage  and  storage  of  the  fluid  (Brown,  1980) .  The  thickness  of  the 
flow  paths  between  the  two  wells  is  believed  to  be  approximately  100 
m,  giving  a  flow  rate  of  about  1  gpm  for  a  unit  depth  of  one  meter. 
This  is  the  test  we  modeled  in  the  application  problems  (Chapter  V) . 

These  tests  demonstrated  that  a  Hot  Dry  Rock  system  can  be 
constructed  and  operated  to  produce  superheated  water  at  temperatures 
suitable  for  generating  electricity.  However,  commercial  power  plants 
require  higher  rates  of  heat  production,  reduced  water  losses  and  a 
credible  basis  for  predicating  useful  lifetime  of  the  heat  source. 


1 . 3  Previous  Work 

Wilson  and  Witherspoon  (1970)  have  reviewed  the  work  on  fluid 
flow  through  fractured  rocks,  in  the  field  of  both  groundwater 
hydrology  and  petrolem  engineering.  Their  work  includes  an  extensive 
set  of  references  that  is  essentially  complete  to  1970.  In  this 
report,  following  the  convention  of  Wilson  and  Witherspoon,  the  term 
"fracture"  is  used  for  most  discontinuities  within  a  rock  mass.  The 
word  "joint"  will  be  generally  employed  in  connection  with  the  finite 
element  joint  model. 

The  general  approach  to  analysis  of  fluid  flow  through  fractures 
has  been  to  model  fluid  flow  through  fractures  assuming  viscous , 
incompressible  flow  between  smooth  parallel  plates  (Snow,  1965).  The 
validity  of  the  cubic  law  for  laminar  flow  of  fluids  through  open 
fractures  consisting  of  parallel  planar  plates  has  been  established  by 
Witherspoon  (1980)  and  Thomas  (1987)  in  laboratory  work.  In  an  "open" 
fracture  the  planar  surfaces  remain  parallel,  and  thus  are  not  in 
contact  at  any  point. 

Deviations  from  the  parallel  plate  model  are  expected  because 
real  joint  surfaces  are  rough  and  contact  each  other  at  discrete 
points.  In  this  thesis,  realistic  rough  surfaces  were  considered 
using  a  factor  of  roughness.  It  has  been  discussed  by  Brown  (1987) 
and  Witherspoon  (1980). 

Witherspoon  and  Noorishad  developed  a  finite  element  model  of 
discrete  fracture  systems.  The  model  coupled  stress  and  fluid  flow 
behavior  in  fractured  rock  masses.  Direct  application  is  to  fluid 
flow  problems  in  hydraulically  fractured  reservoirs  and  naturally 
fractured  rocks  (Noorishad,  Ayatollahi  and  Witherspoon,  1982) . 


A  two  dimensional  finite  element  model  of  fluid  flow  in  fractured 
rock  masses  was  developed  by  Hilber  and  Taylor  (1979).  The 
discontinuities  are  deformable  and  constitute  the  flow  paths.  The 
model  includes  interaction  between  the  fluid  and  the  fracture  motions 
as  well  as  inertia  effects.  They  developed  a  computer  code  based  on 
this  theory.  It  determines  the  hydrodynamic  state  of  the  fluid,  the 
displacement,  strain  and  stress  response  histories  of  the  rock  masses, 
the  change  of  the  kinetic  and  the  potential  energy  of  the  rock,  and 
the  amount  of  energy  dissipated  during  slip. 

A  continuum  approach  has  been  developed  for  modeling  mass 
transport  in  fractured  rocks  by  Schwartz'  and  Smith  (1988).  It 
involves  a  new  application  of  a  particle  tracking  method  in  which 
physical  transport  is  simulated  in  terms  of  velocity  and  the 
variations  in  velocity.  This  model  successfully  duplicates  patterns 
of  anisotropic  dispersion  predicated  by  de  Josselin  de  Jong  (1972). 
Applications  demonstrate  that  the  de  Josselin  de  Jong  approach  for 
estimating  dispersivities  for  idealized  networks  cannot  generally  be 
applied  to  networks  formed  from  sets  of  finite,  irregularly- spaced 
fractures.   This  model  does  not  include  rock  deformation  effects. 

1.4  Objectives 

Our  goal  is  to  develop  a  finite  element  fluid  model  to  simulate 
flow  through  fractured  rock.  This  model  is  the  first  step  in 
developing  a  completely  coupled  flow/rock  deformation/heat  transfer 
model  of  the  Hot  Dry  Rock  reservoir.  Figure  1.2  shows  a  schematic  of 
the  Hot  Dry  Rock  reservoir.  Discrete  flow  paths  are  used  for  the 
analysis,   forming  a   fracture  network  connecting  the  two  wells.   We 


model  a  plane  section  of  the  network,  consisting  of  flow  paths  and 
blocks  of  rock  masses.  Two  types  of  joints  are  present,  shear  joints 
that  are  closed  and  tension  joints  that  are  initially  open.  Fluid  can 
be  stored  in  the  open  joints. 

We  recognize  that  the  density  and  viscosity  of  the  fluid  are 
functions  of  temperature.  In  the  present  model,  we  assume  they  are 
constant,  as  assumed  by  Hilber  and  Taylor  (1979)  in  their  model.  We 
also  assume  the  rock  is  not  deformable  and  that  the  joint  openings  are 
constant.  In  the  future  development  of  this  model,  the  coupling 
between  rock  deformation  and  heat  transfer  will  be  included. 

The  fluid  model  extends  previous  work  by  using  a  solution  scheme 
applicable  to  nonlinear  problems  with  a  large  number  of  unrestrained 
rock  masses.  Another  important  feature  is  the  capability  to  model 
fluid  storage  in  open  joints.  These  joints  aire  filled  as  fluid  is 
pumped  into  the  reservoir.  The  model  includes  the  capability  to 
simulate  the  use  of  tracers  where  tracers  are  pumped  into  the  input 
well  and  the  time  history  at  the  output  is  monitored. 

In  this  thesis,  we  discuss  only  the  fluid  model.  However,  the 
fluid  model  was  developed  as  part  of  a  larger  structural  model,  with 
the  intent  of  coupling  the  two  models.  After  the  coupling,  these 
features  can  give  the  engineer  a  tool  to  simulate  flow  and  to  check 
the  results  using  interactive  computer  graphics.  The  engineer  can 
observe  flow  paths  develop  as  injection  is  started  and  will  be  able  to 
follow  the  flow  until  it  exits  from  the  production  wellbore.  This 
tool  will  aid  in  understanding  and  predicting  the  fluid  behavior  of 
the  Hot  Dry  Rock  reservoir. 
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Figure  1.1:  Schematic  drawing  of  Phase  II  HDR  Reservoir 
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Figure   1.2:    Schematic   drawing  of  HDR  Reservoir  simulation 


CHAPTER  II 
FLUID- FLOW  NETWORK  MODEL 

In  this  chapter,  we  develop  a  model  for  fluid  flow  through  non- 
porous  rock  masses.  The  fluid  flow  model  is  similar  to  the  one 
discussed  by  Hilber  and  Taylor  (1979) .  First  we  derive  the  relation 
between  pressure  gradient  and  flow  rate  for  flow  between  parallel 
plates  (cubic  law) . 

2.1  Derivation  of  flow  model 

We  first  consider  the  fully-developed  laminar  flow  between 
infinite  parallel  plates.  The  plates  are  separated  by  a  distance  a, 
as  shown  in  Figure  2.1.  The  plates  are  considered  infinite  in  the  Z 
direction,  without  variation  of  any  fluid  parameter  in  this  direction. 
The  flow  is  further  assumed  to  be  steady  and  incompressible. 

Because  of  the  non-slip  condition  at  the  wall,  the  X  component  of 
velocity  is  zero  at  both  the  upper  and  lower  plates.  These  two 
boundary  conditions  are  as  follow: 

at   y  =  +a/2        u=0  (2.1) 

y  =  -a/2         u=0. 

Since  the  flow  is  fully  developed,  the  velocity  can  not  vary  with 
X  and  is  a  function  of  Y  only,  u  =  u  (y) .  There  is  no  component  of 
velocity  in  either  Y  or  Z  directions,  v  -  w  =  0. 

A  differential  control  volume  of  size  dV  =  dXdYdZ  is  selected 
(Fig. 2.1)  and  applied  to  the  X  component  of  the  momentum  equation: 


_  J    updV  +   J 


F   +F   = J    updV  +   J    upv.dA  (2.2) 

Sx    Bx   dt     cv         cs 
(1)    (2)      (3)  (4) 

F_   represents  surface  force  and  F_   represents  body  force.   Term 

(4)   indicates   the  momentum  flux  through  the  control  surface  area  dA. 

Term   (3)  represents  the  rate  of  change  of  momentum.   For  steady  flow, 

all  fluid  parameters  are  independent  of  time  and  term  (3)  equals  zero. 

The  net  momentum  flux  through  the  control  surface  is  zero  as  a  result 

of  fully-developed  flow  and  term  (4)  equals  zero.   We  assume  there  is 

no  body  force  in  X  direction,  thus  term  (2)  equals  zero.   Finally,  the 

momentum  equation  reduces  to, 

F   -  0  .  (2.3) 

sx 

We  now  sum  the   forces  acting  on  the  control  volume  in  the  X 

direction.    There  are  normal  forces  (pressure  forces)  acting  on  the 

left  and  right   faces;   there   are   tangential  forces  (shear  forces) 

acting  on  the  top  and  bottom  faces. 

Suppose  the  pressure  at  the  left  face  of  the  element  is  p,  then 

the  force  on  the  left  face  is 

p  dydz  ,  (2.4) 

and  the  force  on  the  right  face  is 

dp 
-  (  p  +     dx)  dydz  .  (2.5) 

dx 


10 


Similarly  for  the  shear  stress   t   ,  the  shear  force  on  the  bottom 

yx 


face  is 


-   r   dxdz  ,  (2.6) 

yx 

and  the  shear  force  on  the  top  face  is 

dr 
yx 

(r   +  dy)  dxdz  .  (2.7) 

yx   dy 

Summing  the   forces  acting  on  each  face  of  the  control  volume,  we 

can  simplify  the  equation  to 


(2.8) 


(2.9) 


dp     yx 

+       =0 
dx    dy 

dr       A 
yx   dp 

—    —  constant. 

dy     dx 

Integrating  this  equation, 

we  obtain 

dp 

r    =      y  +  C.  . 
yx   J          1 

(2.10) 
dx 

Using  the  definition  of  Newtonian  fluid  r    =  /i  (du/dy)  where 

yx 

/i  is  the  dynamic  viscosity,  gives 

du     dp 

m  _  = y  +  c,  (2.11) 


dy     dx 
and 


1   dp   2     1 
u  =       r  v  + 


C 

y  +  C9  .      (2.12) 


2/i   dx 
Applying  boundary  conditions,  gives 
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cx  -  o, 


(2.13) 


1     dp  a 

]2  "   " ' 

H     dx     8 


(2.14) 


and 


1     dp      ,2a. 
u  -  _1_   (y   -   ) 


(2.15) 


2\i     dx 


We  calculate  the  average  velocity  u 


and  obtain 


[a/2udy=   [*/2-  j_  5: 

J_a/2  J -a/2   2/i     dx 


(y2    -  a     )   dy   ,      (2.16) 


a  dp 

u  -  -(  )  _ 

12/i  dx 

The  volumetric  flowrate  is  given  by 


Q  - 


V.  dA 


For  a  unit  depth  1  in  the  Z  direction 

Q 


a/2  - 

ul  dy    , 
-a/2 


or 


(2.17) 


(2.18) 


(2.19) 


Q 
T 


a/2      1        dp  2      af- 

(_)(   y    -   _   )    dy 

_a/2   2\i     dx  4 


(2.20) 


Thus  the  volumetric  flowrate  per  unit  depth  , (q)  is 


a    dp 

q  -  -  ( )  (_)  • 

12/i   dx 


(2.21) 
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In  real  applications,  especially  for  those  small  apertures, 
fracture  walls  are  not  strictly  smooth  and  surface  roughness  can  cause 
turbulence  or  a  boundary  layer  in  the  flow  (Ryan  and  Kimbrell,  1987). 
To  approximate  this,  we  include  a  factor  of  roughness  which  can  cause 
an  apparent  reduction  in  flow, 

a     dp 

q  -  -  ( )  (_)   •  (2.22) 

12/if   3x 

According  to  the  investigations  of  Witherspoon  (1980)  and 
Sundaram  (1987),  the  factor  f  varies  from  1.04  to  1.78  which  depends 
on  the  mechanical  properties  of  fractures  and  rock.  According  to  the 
investigations  of  Louis  (1969)  and  Wilson  (1970) ,  turbulence  only 
exists  in  a  small  portion  of  the  network  which  has  large  hydraulic 
gradients.  Louis  (1969)  states  the  effect  of  turbulence  is  only 
small  within  the  fluid  system,  and  that  the  assumption  of  laminar  flow 
can  be  considered  a  good  approximation.  Only  when  large  portion  of 
the  fluid  network  is  turbulent,  the  total  flow  will  be  significantly 
over  estimated  by  the  laminar  flow.  The  validity  of  the  cubic  law  for 
laminar  flow  through  open  fractures  consisting  of  parallel  plates  has 
been  estabilished  with  apertures  ranging  down  to  a  minimum  of  0.2  um 
(Witherspoon  et  al.,  1980). 

3 
Defining  Kp  =  a  /12juf  as  the  joint  permeability,  we  obtain  the 

relation  between  flowrate  and  pressure  gradient, 


dp 
q  =  -  Kp   (_  )  .  (2.23) 

dx 
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2.2  One  Dimensional  finite  element  model 
2.2.1  The  differential  fluid  element 

We  now  derive  the  differential  equation  for  one  dimensional  fluid 
flow.  The  differential  fluid  element  with  a  unit  depth  is  shown  in 
Figure  2.3.  The  density  of  the  fluid  is  assumed  to  be  constant.  The 
size  of  joint  opening  will  change  the  volume  of  the  differential 
element,   so   it  must  be  included  in  the  formulation.   It  is  assumed 

that  p   is  the  fluid  density,  u  is  the  average  velocity  of  fluid  and  a 
is  the  joint  opening.   From  the  law  of  conservation  of  mass: 

a         a 

pua  -  [pua  +  (pua)dx  J  -  (pa)dx  ,   (2.24) 

ax         at 


we  obtain 


(/ma)  -     (pa)  .  (2.25) 


Since  p   =  constant, 


ax       at 


a      aa 

(ua)  -     .  (2.26) 


ax      at 

From  the  definition  of  volumetric  flowrate  Q  -  ual  (1  is  depth) 
define  the  volumetric  flowrate  per  unit  depth  as  q  =  ua 

aq 

+  a  =  0  ,  (2.27) 

ax 

where  the  derivative  with  respect  to  time  is  indicated  by  ('). 
Assuming   the   flow  rate   is  proportional   to   the  pressure  gradient 
(equation  2.23),  where  K  is  the  joint  permeability, 

we  obtain; 
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3       dp 

(K  )  -  a  -  0  .  (2.28) 


3x       dx 

2.2.2  Quadratic  one -dimensional  element  formulation 

The    interpolation  equation  for  the  one -dimensional   quadratic 
element  (Figure  2.2)  is 

*(£)  -  ax  +  a2C  +  a3|2   .  (2.29) 

The   shape   functions   for  the  one -dimensional  quadratic  element 
relative  to  the  £  -  coordinate  system  are 

Nx  CO  -(£/2)  (£-1) 

N2  (|)  -  -(?-l)  (C+D  (2.30) 


N3  (£)  -(1/2)  (£+1) 


so 


X  (?)  -  Nt(0  Xx  +  N2(|)  X2  +  N3(0  X3  .      (2.31) 


The  Jacobian  of  the  transformation  is 


Y   -  Y 
dx(|)  A3    1 

-  6(3^  -  2X2  +  X3)  +  .     (2.32) 


d|  2 

Assuming  node  2  is  at  the  middle  point  of  the  element, 

Xx  -  2X2  +  X3  -  0  (2.33) 


and 


dX(?)    L 

d  t        2 


(2.34) 
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where  X-  -  X,  -  L,  the  element  length. 

The  row  vector  [D]  contains  the  derivatives  of  element  shape 
functions  with  respect  to  X  written  in  terms  of  the  £  -  natural 
coordinate  system, 

dNx     dN2      dN3 

[D]  -  {  (O   (O   (O  )  .         (2.35) 

dX      dX       dX 

We  use  the  chain  rule  to  obtain, 


dNx    dNx  dX 


(2.36) 


d|     dX   d£ 
giving 


2    r         l  l   \ 

[D]  -  {(£  -  )   -2£    (£  +  )}  .    (2.37) 

L   L      2  2  J 


2.2.3  Galerkin's  method 

The  weighted  residual  integral  is  the  line  integral  along  the  X 
direction.  We  will  assume  the  weighting  function  for  the  Sth  node, 
W  ,   consists   of  the   shape   functions   associated  with  the  Sth  node 

(Segerlind,   1984).    Therefore   [W  ]   =   [N  ] ,   [N]  is  the  row  vector 

containing  the  element  shape  functions.   Defining  a  column  vector  {R} , 
each  component  of  (R)  represents  a  residual  equation, 


(R) 


T    3                 dp 
[N]  {  (Kp  )  -  a  }  dx  ,      (2.38) 


L       3x       dx 
separating  the  terms  inside  the  integral 
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<R) 


t  a     ap 

[N]  (Kp  )  dx 


ax 


ax 


[N]  a  dx  .      (2.39) 


Integrating  by  parts,  gives 


T   3p 
{R}  =  [N]  Kp  

ax 


dp  3[N] 

Kp dx 

L    dx   3x 


[N]  adx  .   (2.40) 


The  integration  over  the  body  in  equation  (2.40)  is  performed 
dividing  the  body  into  elements  and  summing  the  integration  over  each 
element.   For  each  element,  we  have, 

3p 


ax 


[D]  [P] 


(2.41) 


We  can  obtain 


{R}  -  [N]  q 


Kp[D][P][D]  dx 


[N]  a  dx  .(2.42) 


Because  a  =  [N]  [a]  and  both  [P]  and  [a]  are  constant  with  respect 
to  the  integration, 

L, 


{R}  =  [N]  q 


Kp[D][D]  dx[P]  - 


[N]  [N]  dx[a]  .  (2.43) 


From  the   element  matrix  formulation  equation  (2.37),  we  expand 

T 
[D]  [D]  to  obtain, 


[D]T[D] 


(£-1/2)* 

-2?(M/2) 


2^-1/2) 


(r_i/4) 

-2£(£+l/2) 


(|  -1/4)    -2£(£+l/2)    (1+1/2)^ 
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[B] 


(2.44) 


We  use   the   shape   functions  to  describe  the  fracture  opening  in 
3 


equation  2.23,  Kp  -  a  /12/if 

-  {[N]  [a]}3/12/if 


(2.45) 


substitute  (2.44)  and  (2.45)  into  equation  2.43,  gives, 

L 

i  n 

3 


{R}    =    [N]    q 


(3) 


L0         12/il/f  . 


im    [a]KfB]dx   [P]    - 


(2) 


[N]  [N]dx  [a](2.46) 


(1) 


Integration  is  done  numerically  by  transforming  to  the  natural 
element  coordinate  system.  We  can  transform  term  (1)  from  the  global 
coordinate  system  X  to  the  natural  coordinate  system  £  by, 


IT. 

[N]  [N]   d£  [a]=  [S]  [a]  , 
-1 


where  [S]  is  the  storage  matrix  due  to  crack  opening  velocities 
Similarily,  transforming  the  stiffness  matrix, 


(2.47) 


3/i  L2f 


{[N]     [a]T    [B]    dx  -  __^_ 
L  6/i   L   f 


1([N]    [a]}3    [B]    d£.(2.48) 

-1 


Using  Che  shape  functions  in  terms  of  the  natural  coordinates, 
we  obtain  for  the  stiffness  matrix, 


18 


[K] 


6/iLf 


1[(al      -    a2  +  ^3__)i2  +    (^   -   a!_)g   +a2]3TB]    d£,(2.49) 
-12  2  2  2 


Finally  we   get, 

{R}  -  [Q]  -  [K][P]  +  [S][a]  .  (2.50) 

Setting  {R}  =  {0} 

[K][P]  =  [Q]  +  [S][a]  ,  (2.51) 

where  [K]  is  the  permeability  (stiffness)  matrix,  [P]  is  the  column 
vector  containing  nodal  pressures,  [Q]  is  a  vector  of  specified  flow 
rates,   [S]   is   the  storage  matrix  due  to  the  opening  of  the  gap,  and 

[a]  is  the  gap  opening  velocity  at  the  nodes.  We  can  solve 
equation  (2.51)  for  nodal  pressures. 

2.3  One  Dimensional  Tracer  Model 

Tracer  are  used  to  track  the  motion  of  fluid  in  the  system  of 
fractures  and  can  provide  important  information  on  the  fractures 
system.  Typically,  the  tracer  is  a  chemical  or  radioactive  material 
that  is  mixed  with  the  injected  fluid  and  the  concentration  is 
monitored  at  the  output.  In  Fenton  Hill  experiments,  engineers  have 
used  sodium  fluorescein  dye,  sodium  bromide  and  sodium  nitrate  to 
estimate  the  volumes  of  fractured  systems. 

The  first  derivation  of  the  tracer  model  is  for  an  active  element 
already  filled  with  fluid.  The  differential  tracer  element  with  a 
unit   depth   is   shown   in  figure   2.4.    The   fluid   is   assumed 

incompressible.   The  concentration  of  tracer  is  denoted  by  c ,  u  is  the 
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average  velocity  of  fluid  and  a  is  the  gap  opening.   From  the  law  of 
conservation  of  mass: 


3  3 

cua  -  [cua  +  (cua)dx  II  -  (ca)dx  .      (2.52) 

3x  at 


Differentiating  the  right  side, 


3  3c        3a 

(cua)   -  a  +  c  .  (2.53) 


ax         at      at 

The  differential  form  becomes 


3c      1  a  3a 

=  (    -     (cua)    -   c    )  .        (2.54) 


3t      a         3x  3t 

We  solve  equation  2.54  using  a  simple  finite  difference  scheme. 
It  is  assumed  that  a  fluid  element  is  divided  into  three  sections . 
For  each  section,  the  current  tracer  is  going  to  next  section  and  will 
be  replaced  by  the  new  coming  tracer.  The  movement  depends  on  the 
direction  of  fluid  velocity.  Equation(2 . 54)  gives  the  rate  of  change 
of  concentration  in  each  section  of  fluid  element. 


c      c 
3c      new    old 


(2.55) 


3t  At 


and 


S  ,  S  CUa   I  •     _   CUa   I  /O   C/-N 

(cua)  =■  J_in 'out  (2.56) 

3x  Ax 
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da  a    -  a  ,  , 

new    old 


(2.57) 


at       At 

At  cua  I .   -  cua  I         a    -  a  ,  , 
'in       'out       new    old 

c    -  c    +  ( -c  )   (2.58) 

new    old    a         Ax  At 


Equation(2 .58)  shows  the  relationship  between  current  tracer  and 
the  new  coming  tracer. 

In  the  fluid  system  network,  mixing  of  the  tracer  occurs  at 
the  connecting  nodes  of  fluid  elements.  The  calculation  of  tracer 
mixing  is  illustrated  in  Figure  2.5,  where  there  are  two  branches  of 
fluid  flowing  into  a  fluid  element.  At  node  I  the  calculation  of 
concentration  of  input  tracer  is  a  weighted  average  of  the  input 
concentration, 


u..  c..  +  u2  c2 

c.     -  .  (2.59) 

input   

ul  +  u2 


The  output  concentration  is  the  current  concentration  at  the  node. 
For  nonactive  fluid  elements,  the  calculation  of  mixing  of  tracer  is 
as  follows  (refer  to  Figure  2.5).  The  input  node  for  any  inactive 
element  may  be  connected  to  several  other  active  or  nonactive 
elements .  The  active  elements  give  the  total  flow  into  the  node .  The 
flow  into  each  inactive  element  is  distributed  using  the 
permeabilities  of  the  inactive  elements. 

In  Figure  2.5,  there  are  three  inactive  elements  at  the  input 
node.   The  input  to  element  number  23  is  then: 
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At  node  I 


and 


K    . 
Pi 

qin=    (    *1  +  ^2    > 

K    .    +  K  .,    +  K  „ 
pi          pi          p2 

K    . 

Pi 

Cin=    (    cl  +  c 2   > 

(2.60) 


(2.61) 


pi    pi    p2 


The  concentration  of  tracer  in  this  element  becomes 


c.  Q.  At  +  c,  Q.  At  +  c  ,  .  V,..,,  , 
k  xk       old  filled 

Q.  At  +  Q,  At  +  V-..,  , 
xi      xk       filled 


c„ld- •        <2-62> 


K    indicates   the  conductivity  of  fluid  element,   q   is   the 

flowrate   in  the  active  element,  and  c  is  the  concentration  of  tracer. 

V.c--,t  ,  means   the  volume  of  the  element  which  is  already  filled  with 
filled  J 

fluid. 
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Figure  2.1:  Fluid  flow  between  two  parallel  plates 
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Figure  2.2:  One -Dimensional  Quadratic  element 
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Figure  2.3:  Differential  fluid  element 
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Figure  2.4:  Differential  tracer  element 
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Figure  2.5:  Example  of  tracer  mixing  at  a  joint 
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CHAPTER  III 
DESCRIPTION  OF  FLUCRK 

This  chapter  provides  a  brief  description  of  FLUCRK.  Significant 
features  of  the  implementation  include  the  use  of  dynamic  relaxation, 
the  connection  between  the  structure  model  and  fluid  model,  the  use  of 
active/nonactive  elements  to  model  storage,  and  the  use  of  interactive 
computer  graphics  to  specify  the  problem  and  examine  the  results.  The 
flow  chart  of  FLUCRK  is  shown  in  Fig. 3.1. 

3.1  User  Interaction 

3.1.1  Automated  mesh  generation 

Before  executing  FLUCRK,  the  user  must  generate  a  mesh.  A  short 
program  called  "AU-MESH"  allows  the  user  to  generate  simple  symmetric 
meshes  interactively.  The  program  automatically  creates  an  output 
geometry  file  which  can  be  read  by  FLUCRK.  This  program  was  used  to 
generate  the  meshes  used  in  application  calculations  (chapter  V) . 

3.1.2  Interactive  Problem  specification 

After   the   geometry  file  has  been  read  into  FLUCRK,  the  problem 

specification   is  performed  interactively.    The  user  can  modify  the 

mesh  by  deleting  or  adding  elements.   The  user  can  enter  the  FLUID  MAT 

COND  page   to   examine  and  change  material  properties   interactively. 

The   FLUID  BC  page   allows   the  user  to  specify  boundary  conditions 

(either  pressure  vs   time  or  velocity  vs  time)  at  arbitrarily  picked 

nodes.    The  far  field  boundary  conditions  (Section  3.4)  are  specified 

Vi 
using  a   flowrate  vs  pressure  curve.   In  the  ANALYSIS  PARAMETER  page 
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the  user  can  set  the  start  time,  the  end  time,  the  time  intervals  to 
store  analysis  results,  restart  data,  and  plot  data.  Parameters 
controlling  the  solution  are  specified  in  the  FLUID  page. 

Besides  these  parameters,  the  user  can  also  specify  the  gap 
opening  used  to  flag  whether  an  element  is  active.  Elements  with 
large  initial  openings  accumulate  fluid  until  they  are  filled, 
elements  with  smaller  gaps  are  assumed  filled  and  initially  active. 

At  the  MAIN  page,  there  are  two  options  used  to  manually  save 
problem  data.  These  are  the  SAVE  RESTART  and  SAVE  PLOT  functions.  SAVE 
RESTART  saves  all  problem  data  at  the  present  time  in  a  restart  file. 
The  user  can  then  retrieve  all  the  problem  data  to  rerun  the 
analysis.   SAVE  PLOT  will  save  plot  data  for  the  current  time. 

Once  a  problem  is  finished,  the  user  can  always  return  to  change 
parameters  or  boundary  conditions  and  execute  the  problem  using  the 
RESTART   function  or  examine  the  results  by  entering  the  PLOTS  page. 

In  the  FLUID  PLOT  page,  there  are  nine  functions  to  examine 
results  already  obtained.  These  functions  are  divided  into  two  parts. 
One  is  color  plots  of  a  specified  parameter  (pressure,  flowrate, 
velocity  or  tracer)  at  a  chosen  time  step.  Another  is  time  history 
plots  (pressure,  flowrate,  velocity,  tracer  or  kinetic  energy). 
Quantities  are  displayed  as  a  function  of  time  at  a  specified  node  or 
element . 
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3.2  Connection  Between  Structure  Elements  and  Fluid  Elements 

This  finite  element  model  was  developed  as  part  of  a  large 
structural  code,  CRACKER  (Swenson,  1985),  which  can  solve  both  dynamic 
and  static  problems.  CRACKER  uses  six  noded  triangular  elements  to 
model  the  continuum  and  includes  interface  elements  to  model  contact 
between  bodies.  We  added  a  third  element  type,  a  one  dimensional 
fluid  element.  As  illustrated  in  Figure  3.3,  each  fluid  element  lies 
between  two  structural  elements  and  has  six  nodes.  However,  only 
three  nodes  per  element  are  required  to  represent  the  network  of  fluid 
elements.  To  save  storage  in  the  fluid  solution,  a  fluid  model 
numbering  system  is  used  that  is  independent  of  the  structural  nodes. 
The  fluid  numbering  system  (P, ,  P« ,  P~,  etc)  is  invisible  to  the  user, 

but  is  automatically  generated  when  the  fluid  solver  is  entered. 

We  create  two  array  listings  to  define  the  connection  between  the 
fluid  and  structure  elements.  The  first  array  is  called  a 
connectivity  list  and  indicates  the  node  numbers  of  the  structure  and 
fluid  nodes  in  each  fluid  element.  We  use  the  following  storage, 
where  ELEM  indicates  the  fluid  element  number, 

C0NN(1-6,ELEM)  :  STRUCTURE  NODE  NUMBERS 
C0NN(7,ELEM)    :  ELEMENT  TYPE 

C0NN(8,ELEM)    :  MATERIAL  NUMBER  (SPECIFIED  BY  USER) 
C0NN(9-11,ELEM):  FLUID  NODE  NO. 
This   array  allows   the  program  to   go  from  the  fluid  elements 
to  the  fluid  nodes . 

Another  array  listing  which  we  call  the  inverse  connectivity  list 
stores   the  number  of  all  fluid  elements  connected  to  each  fluid  node. 
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This  array  gives  the  program  pointers  from  the  fluid  nodes  to 
elements.  The  array  is  used  when  calculating  net  flow  into  any  node. 
For  a  fluid  node,  FLUID,  this  data  is  stored  as, 


FL0W_C0N(1, FLUID) 
FL0W_C0N(2,FLUID) 
FLOW  PT( FLUID) 


NO.  OF  FLUID  ELEMENTS  CONNECT  TO  THE  NODE 
POINTER  TO  ELEMENT  NUMBER  LIST 
ELEMENT  NUMBER  LIST 
The  above  data  storage  scheme  eliminates  the  need  for  repeated 
searching  to  identify  connectivity. 

3.3  Active/Nonactive  Element 

Conceptually,  the  joint  network  includes  joints  that  are  initially 
open.  These  joints  must  be  filled  with  fluid  before  we  can  assume 
that  conservation  of  mass  applies  for  that  joint.  Until  the  open 
joints  are  filled,  they  are  called  nonactive.  The  pressure  at  the 
nodes  of  nonactive  elements  are  assumed  zero  and  they  are  not  included 
in  the  solution  of  equation  (2.51).  Instead,  they  supply  boundary 
conditions  to  this  equation.  Once  it  is  full,  the  element  is  turned 
active  and  included  into  the  pressure  calculation. 

For  all  initially  inactive  elements,  the  initial  empty  volume 
is  calculated.  In  each  time  step,  the  flow  into  the  inactive  elements 
is  calculated  and  the  new  empty  volume  obtained.  The  total  flow  into 
the  empty  joint  is  calculated  and  used  to  obtain  the  new  empty  volume, 


vnew  m   vold   _  (Total  flow  rate  in)  *  TSTEp 
empty    empty 
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where  TSTEP  is  the  time  increment.   If  V      is  smaller  than  or  equal 

empty  M 

to   zero,   that  means   the  fluid  element  is  filled  with  fluid.   It  is 

automatically  turned  on  to  active  status  and  included  into  the  flow 

calculation. 

In  general,  several  active /no nactive  elements  may  connect  at  one 
node.  Total  flow  into  the  node  is  obtained  using  all  active  elements. 
If  more  than  one  nonactive  element  is  connected  to  a  node,  the  flow  is 
proportioned  to  the  nonactive  elements  using  their  conductivities. 

If  an  element  fills  in  the  middle  of  a  time  step,  all  additional 
flow  into  that  element  is  lost  to  the  calculation.  It  is  important  to 
choose  the  time  step  small  enough  that  this  loss  is  minimized. 

3.4  Far  Field  Boundary  Condition 

The  finite  element  model  can  not  extend  to  infinity,  including  all 
flow  paths  between  the  wells.  Instead,  we  can  include  only  a  finite 
region  around  the  two  wells.  We  call  this  boundary  around  the  problem 
the  far  field  boundary.  In  the  real  problem,  flow  leaks  out  of  this 
boundary  into  the  infinite  system.  As  a  first  approximation,  we 
assume  the  leakage  is  a  function  of  pressure  at  the  far  field  boundary 
nodes.  The  user  can  specify  a  nonlinear  relation  of  pressure  vs 
f lowrate .  As  will  be  discussed  in  chapter  V,  we  examine  the 
sensitivity  of  the  solutions  to  the  far  field  boundary  conditions. 

3.5  Solution  Technique 

Dynamic  relaxation  (Day,  1965;  Underwood,  1983;  Papadrakakis , 
1981)   is   the  technique  used  to  obtain  the  solution  to  equation  2.51. 
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Dynamic  relaxation  (DR)  is  an  explicit  iterative  method  for  the  static 
solution  of  simultaneous  equations.  It  is  based  on  the  fact  that  a 
system  undergoing  damped  vibration  excited  by  a  constant  force 
ultimately  comes  to  rest  in  the  displaced  position  of  static 
equilibrium  of  the  system.  It  has  been  extensively  used  for  both 
linear  and  nonlinear  structure  analysis  (Frieze,  Hobbs  and 
Dowling, 1978 ;  Pica  and  Hinton,  1980).  This  method  is  derived  from  the 
dynamic  equilibrium  equation: 

[M][P]   +[C][P]  +[K][P]  -  [Q]  ,         (3.1) 

where  [M]  is  a  mass  matrix,  [C]  is  a  damping  matrix,  [K]  is  the 
stiffness  matrix,  [P]  is  the  fluid  pressure  vector,  [Q]  is  the 
flowrate  vector  and  the  dot  represents  differentiation  with  respect  to 
time. 


Using  the  method  of  central  difference,  we  can  get  [P] ,  [P]  and  [P] 
as  follows  : 


. .   [QJ-CP-KP 
P  , 

M 


P    =  P    +  P  h  ,  (3.2) 

new     old 


P    -  P,.  +  Ph 
new     old 


where  h  is  a  fixed  time  increment.  To  preserve  the  explicit  form  of 
the  central  difference  integrator  [M]  must  be  diagonal  and  the  damping 
matrix  [C]  has  the  form  [C]  -  c  [M] . 


31 


To  obtain  the  static  solution  from  the  equation,  we  should  select 
the  damping  coefficient  c,  the  time  increment  h,  and  the  mass  matrix 
[M]  to  converge  rapidly.  Physically,  the  damping,  time  increment  and 
mass  are  selected  so  that  the  transient  response  will  reach  the  steady 
state  under  the  applied  load  Q.  Only  P  and  Q  must  represent  the 
physical  problem  and  c  and  M  do  not  need  to  represent  the  physical 
structure  (Underwood  1983) . 

The  pseudo  mass  matrix  [M]  here  being  used  is  derived  from 
Gerschgorin' s  theorem  : 

1    2 

(3.3) 

where  m. .  are  the  diagonal  elements  of  the  pseudo  mass  matrix  and  K. . 

are  the  elements  of  the  stiffness  matrix  K.  Underwood  (Underwood 
1983)  suggested  evaluating  [M]  for  h  -  1.1  and  iterating  with  h  =  1 
to  ensure  stability.  We  found  it  was  necessary  to  use  h  =  0.5  for 
stability. 

The  damping  matrix  coefficient  is  computed  at  each  iteration  from 
Rayleigh's  quotient  as  (Underwood,  1983)  : 


1 

2 

m       -   — — 
ii          4 

• h  z 

K 

■/ 


P]T[K]  [P]  /  [P]T[M]  [P]  (3.4) 


The  flow  chart  of  dynamic  relaxation  is  shown  in  Figure  3.2. 

We  have  tested  two  different  convergence  criteria.  The  first  uses 
the  maximum  kinetic  energy  during  the  solution.  At  each  iteration, 
the  current  kinetic  energy  is  compared  to  the  maximum  kinetic  energy 
times  a  factor.  If  the  current  energy  is  smaller,  we  assume  that  the 
system   is   converged.    The   difficulty  with  this   criteria  is  that 
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pressure  boundary  conditions  can  lead  to  a  large  initial  kinetic 
energy,  and  convergence  is  not  obtained  unless  a  very  small  factor  is 
used.  Under  this  condition,  the  result  is  not  reliable.  The  second 
criteria  checks  the  difference  between  the  previous  pressures  and 
current  pressures  divided  by  previous  pressures.  If  the  difference  is 
smaller  than  the  tolerance  specified  by  the  user,  we  assume  the 
pressures  are  converged.  The  reasonable  tolerance  is  picked  by  a 
series  of  trial  and  error  tests. 

The  main  advantage  of  dynamic  relaxation  lies  in  its  basic 
simplicity  and  the  straightforward  algorithms  which  can  be  written  for 
complicated  problems.  In  the  implementation,  no  global  matrices  are 
actually  formed.  Instead,  equation  3.1  is  solved  by  looping  in  all 
elements  and  accumulating  the  unbalanced  forces.  There  is  no  required 
numbering  scheme  for  the  elements  or  nodes,  which  is  ideal  for  the 
present  application,  that  includes  active  and  inactive  elements  that 
change  during  the  analysis.  An  additional  advantage  is  the 
exceptional  robustness  of  this  solution  scheme  for  nonlinear  problems, 
such  as  when  the  structural  model  is  coupled  to  the  fluid  model  and 
fluid  flow  becomes  a  cubic  function  of  joint  opening  (equation  2.22). 

The  main  disadvantage  of  this  method  is  the  relatively  long 
solution  time  compared  to  direct  solution  methods  for  a  linear  system 
of  equations. 
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Figure   3.2:    Flow  chart   of  DR 
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CHAPTER  IV 


VERIFICATION  PROBLEMS 


A  number  of  verification  problems  have  been  solved  to  test  the 
computer  program  and  to  demonstrate  the  method  of  analysis  described 
in  chapter  II.  Problems  one  and  two  are  designed  to  test  the  effect 
of  fluid  flow  while  all  fluid  elements  are  active  (problem  one)  and 
nonactive  (problem  two)  at  the  beginning.  They  use  the  same  mesh  with 
a  straight  flow  path  that  includes  25  fluid  elements.  Tracer  was 
pumped  into  the  flow  path  in  these  two  problems  to  verify  the  motion 
of  tracer  through  time.  Problem  three  demonstrates  the  mixing  of 
fluid  flow  and  the  change  of  concentration  of  tracer  in  multiple  flow 
paths.  Problem  four  tests  the  specification  of  nonlinear  far  field 
boundary  conditions . 

4.1  Problem  One  -  Flow  in  Crack 

This  problem  is  a  simple  one  dimensional  problem  (Figure  4.1). 

Fluid  flow  was  pumped  into  the  left  side  under  constant  pressure  and 

output   from  the   exit   in  right   side.   Twenty  five  fluid  elements 

construct  a  straight  flow  path  and  they  are  all  filled  with  fluid  at 

the  beginning.    The  problem  was   solved  for  25  time  steps  of  1  sec 

each. 

-5     2 
All   elements  had  the  same  dynamic  viscosity  u  =  0.15x10  N-s/m 

and  gap   opening  a  =  0.03  m,   as   shown  in  Figure  4.1.   The  input 

pressure  was   20   Pa  and  output  pressure  was  0  Pa.   Both  were  kept 
constant  during  the  analysis. 
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The  calculated  flowrates  at  the  input  and  output  elements  are 
shown  in  Figure  4.2  and  Figure  4.3.  They  both  are  identical  in  every 
time  step  and  correspond  to  the  expected  values. 

To  demonstrate  calculations  using  the  tracer  option,  we  injected 
tracer  as  a  pulse  from  2  through  5  sec.  The  input  concentration  was 
100%.  Time  history  plots  of  the  tracer  at  the  input,  center,  and 
output  nodes  are  shown  in  Figure  4.4.  The  flow  of  the  tracer  from  the 
input  to  output  is  clearly  seen,  with  the  pulse  travelling  at  the 
velocity  of  the  fluid.  Some  numerical  diffusion  of  the  tracer  is  seen 
to  occur. 

This  example  shows  the  pressure  decreased  linearly  from  the  input 
to  output  pressure,  and  flowrates  and  velocities  are  constant  through 
the  problem.  It  also  demonstrates  that  the  tracer  calculation  works 
correctly. 

4.2  Problem  Two  -  Flow  in  Empty  Crack 

Problem  two  uses  the  same  geometry  and  mesh  as  problem  one.   The 

only  differences   are  that  all  fluid  elements  (except  those  connected 

to   the   input  and  output  nodes)  are  initially  nonactive  and  the  gap 

opening   is   0.05  m.    This   example   demonstrates   filling  of  empty 

elements . 

2 
The   specified  input  pressure  was  2  N/m  and  the  output  pressure 

2 
was   0  N/m  .    Flowrates   of  the  input  element  and  output  element  are 

displayed  in  Figure  4.5  and  Figure  4.6.   Initial  flowrates  are  high 

due   to   the   empty  elements  next  to  the  input  elements.   As  the  flow 

progressively  fills  the  elements,  the  resistance  to  flow  increases  and 

the   flow  rate   drops.    Figure  4.7  demonstrates   the  pressure  and 
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flowrates  at  middle  of  the  flow  path.  The  middle  element  turned  on 
after  fluid  flow  filled  it  after  9  time  steps.  After  twenty  two  time 
steps,  all  fluid  elements  were  filled  with  fluid,  pressures  decreased 
linearly  between  input  and  output,  and  flowrates  were  constant 
everywhere. 

The  expected  flowrate  is: 


a  dp         ,  (4.1) 


12/xf     dx 


3 
0. 14  m  /  sec 


which  is  identical  to  Figure  4.6. 

Figure  4.8   displays   the  motion  of  tracer  at  input,  middle  and 
output  locations  of  this  crack  system. 

4.3  Problem  Three  -  Flow  in  Simple  Network 

Problem  three  has  ten  fluid  elements  initially  active  with  a  gap 

opening  of  0.01  m  as  shown  in  Figure  4.9.   Fluid  was  pumped  into  two 

2 
entries   in  the   left  side  at  5  N/m  and  output  from  two  exits  in  the 

2 
right  side  at  0  N/m  .   The  calculated  flowrates  at  the  entry  and  exit 

are   displayed  in  Figure  4.10  and  Figure  4.11.   Figure  4.12  shows  the 

fluid  properties  at  middle  flow  path  where  the  pressure  equals  half  of 

the   input  pressure   and  the   flowrate  equals   the  sum  of  two  input 

flowrates . 

The  motion  of  tracer  is  shown  in  Figure  4.13.   Tracer  was  pumped 

into   one   entry  at  a  concentration  of  100%.   The  concentration  of 
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tracer  at  the  middle  of  flow  path  is  only  50%  because  of  the  mixing  of 
the  two  input  flows.  When  the  flow  splits  at  the  output,  the 
concentration  remains  at  50%.  Because  of  the  relatively  coarse  mesh 
used,  diffusion  of  the  tracer  occurs  at  the  output.  Figure  4.14  shows 
the  pressure  plot.  The  pressure  decreased  linearly  from  the  input  to 
the  output  within  each  element.  The  overall  tracer  plot  is  displayed 
in  Figure  4.15,  which  shows  the  reduction  in  concentration  of  tracer 
after  mixing. 

4.4  Problem  Four  -  Far  Field  Boundary  Conditions 

This  problem  tests  the  option  to  specify  a  pressure- flowrate 
relation  at  the  far  field  boundary  conditions  (Figure  4.16).  It  is 
assumed  the  far  field  condition  is  a  nonlinear  relation  between 
pressure  and  flowrate  as  shown  in  Figure  4.17.  The  input  flowrate 
boundary  condition  at  fluid  node  5  is  between  fluid  element  2  and  3. 

3 

The  value  of  the  input  flowrate  is  1  m  /sec. 

Figure  4.18  and  Figure  4.19  show  the  results  of  flowrates  of  fluid 
elements   2   and  3,   with  the  absolute  sum  of  them  equal  to  the  input 

3 
flowrate   1  m  /sec.    Pressures   at  both  ends   and  input  node  are 

displayed  in  Figure  4.20. 
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Figure  4.1:  Verification  problem  one 
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Figure  4.2:  Flowrate  of  input  element  (problem  one) 


X=   98  .  868  Y=     9  . 659 


s  a  see 


e 

see 

LU 

a 

see 

f- 

<X 

C£ 

Z3 

a 

i  ee 

o 

e  eee     e  see     1  eee     i  see     2  eee 

TIME 


Figure  4.3:  Flowrate  of  output  element  (problem  one) 
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Figure  4.4(a):  Tracer  at  input  node  (problem  one) 
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Figure  4.4(b):  Tracer  at  middle  of  crack  (problem  one) 
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Figure   4.4(c):    Tracer  at   output  node    (problem  one) 
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Figure  4.5:  Flowrate  of  input  element  (problem  two) 
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Figure  4.6:  Flowrate  of  output  element  (problem  two) 
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Figure  4.7(a):  Pressure  at  middle  of  crack  (problem  two) 
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Figure  4.7(b):  Flowrate  at  middle  of  crack  (problem  two) 
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Figure  4.8(a):  Tracer  at  input  node  (problem  two) 
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Figure  4.8(b):  Tracer  at  middle  of  crack  (problem  two) 
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Figure  4.8(c):  Tracer  at  output  node  (problem  two) 
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Figure  4.10:  Flowrate  of  input  element  (problem  three) 
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Figure  4.11:  Flowrate  of  output  element  (problem  three) 
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Figure  4.12(a):  Pressure  at  middle  of  crack  (problem  three) 


X=    3  .  897  Y=     4  .  966 


«  9     600  - 


LU 


a  see     i  see     2  aee     3  eae     «  aea 


TIME 


Figure  4.12(b):  Flowrate  at  middle  of  crack  (problem  three) 
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Figure  4.12(c):  Velocity  at  middle  of  crack  (problem  three) 
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Figure  4.13(a):  Tracer  at  input  node  (problem  three) 


51 


X  =  5.009       Y =  5.000 


IXI 

o 

<I     a    298-- 

S- 


e   988  i    aee  2   eee 


TIME 


-t- 


H 


3    see  4    ee>i 

i 

< « 10        > 


Figure  4.13(b):  Tracer  at  middle  of  crack  (problem  three) 
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Figure  4.13(c):  Tracer  at  output  node  (problem  three) 
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Figure  4.14:  Pressure  plot  of  problem  three 


Figure  4.15:  Tracer  plot  of  problem  three 
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Figure  4.16:  Verification  problem  four 
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Figure  4.17:"  Far  field  boundary  condition 
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Figure  4.18:  Flowrate  at  left  end  of  fluid  element  (problem  four) 
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Figure  4.19:  Flowrate  at  right  end  of  fluid  element  (problem  four) 
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Figure  4.20(a):  Pressure  at  left  end  (problem  four) 
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Figure  4.20(b):  Pressure  at  input  node  (problem  four)) 
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Figure  4.20(c):  Pressure  at  right  end  (problem  four) 
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CHAPTER  V 
APPLICATIONS 

We  analyzed  five  problems  which  simulate  the  Phase  II  HDR  project 
of  Los  Alamos  national  laboratory.  Problems  one,  two,  three  and  four 
all  include  630  fluid  elements  and  3339  structural  nodes.  Problem 
five  has  165  fluid  elements  and  924  solid  nodes.  Problems  one,  two 
and  three  compare  the  effect  of  different  far  field  boundary 
conditions.  All  elements  are  active  at  the  beginning,  and  there  is  no 
flow  leakage  (problem  one)  and  zero  pressures  (problem  two)  in  far 
field.  In  problem  three,  a  nonlinear  pressure/flowrate  relation  is 
specified  at  the  far  field  boundary.  In  problem  four,  we  examine  the 
effect  of  an  unsymmetric  flow  path  with  increased  conductivity.  The 
final  problem  displays  the  fluid  motion  in  which  half  of  the  fluid 
elements  are  nonactive  initially.  The  far  field  conditions  of 
problems  three,  four  and  five  are  all  the  same,  which  is  shown  in 
Figure  5.14.   The  input  pressure  and  output  pressure  are  assumed  to  be 

2.034xl07  N/m2  (2950  psi)  and  0  N/m2  (0  psi)  for  all  the  five 
application  problems,  giving  the  same  pressure  difference  as  observed 
experimentally.  Fixing  the  output  pressure  at  zero  helps  speed 
convergence . 

The  dimensions  of  Phase  II  reservoir  are  about  400  m  long  and  200 
m  wide.  It  is  believed,  that  in  a  plan  view,  the  joints  are  spaced  at 
about  10  m  horizontally  and  10  m  vertically.  To  speed  the  problem 
solution,  we  used  a  spacing  of  20  m  between  joints.  The  large  joints 
are  assumed  to  have  an  opening  of  0.004  m  and  the  vertical  joints  have 
an  opening  of  0.0000396  m  (Figure  5.2).   These  values  were  chosen  to 
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give  what  we  believed  to  be  approximately  the  correct  storage  capacity 
and  approximately   the  observed  f lowrate .   The  depth  of  the  flow  path 

-3 
between  wells  is  about  100  m  and  the  total  f lowrate  is  about  6.3x10 

3  -5 

m  /sec   (100   gpm) .    For   one   meter   depth,  the  flowrate  is  6.3x10 

3 
m  /sec  (1  gpm).   From  the  steam  tables  (Haar  and  Gallagher,  1984),  the 

-4      2  o 

dynamic   viscosity  of  water  is  1.416x10   N-s/m  at  350  bar  and  200  C. 

In  our  simulation,  we  included  the  factor  of  roughness  -  1.5  into  the 

stiffness  calculation. 


5.1  Problem  One  -  Zero  Flow  in  Far  Field 

In  problem  one  we  assume  that  there  is  no  leakage  in  far  field. 

The  input  flowrate  and  output  flowrate  are  both  shown  in  Fig. 5. 3. 

After   7500   iterations,  the  system  reached  a  stable  condition  and  the 

-5   3 
output  flowrate  equaled  input  flowrate  at  6.3x10   m  /sec  (1  gpm). 

Figure   5.4   shows   the   pressure  at  a  far  field  node  which  has  no 

flow  leakage.    Figure   5.5   demonstrates   the  convergence  of  kinetic 

energy.      Figure   5.6   shows  the  pressure  distribution.   The  maximum 

pressure  is  represented  by  red  and  the  minimum  pressure  is  represented 

by  green.   The  pressure  decreased  from  the  input  node  to  output  node. 

Figure  5.7  shows  the  flowrate  plot.   The  arrow  indicates  the  direction 

of   fluid  flow.   The  flowrates  of  the  far  field  nodes  are  smaller  than 

1.0x10     which   approximately  equals  to  zero  because  of  no  far  field 
leakage . 
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5.2  Problem  Two  -  Zero  Pressure  in  Far  Field 

The  only  difference  between  problem  one  and  two  is  the  far  field 
condition.  The  specified  pressure  vs  flowrate  far  field  boundary 
condition  was  chosen  to  approximate  a  zero  pressure  condition.  This 
is  shown  in  Figure  5.8. 

Figure  5.9  shows  the  flow  at  the  input  and  output  wells.  The 
input  flowrate  is  approximately  three  times  the  flowrate  of  problem 
one  due  to  reduced  resistance,  because  of  the  added  flow  paths  to  the 
far  field.  The  output  flowrate  is  only  1/4 th  the  output  flowrate  of 
problem  one,  since,  due  to  far  field  leakage,  most  flow  does  not  reach 
the  output  node . 

The  pressure  and  the  leakage  flowrate  at  a  far  field  node  are 
shown  in  Figure  5.10.  The  pressure  is  approximately  zero.  Figure 
5.11  shows  the  convergence  of  kinetic  energy. 

Figure  5.12  shows  the  pressure  distribution.  Far  field  pressures 
are  approximately  zero.  The  flowrate  plot  is  shown  in  Figure  5.13. 
The  arrows  of  flowrates  in  far  field  elements  point  outward  which 
indicates  flow  leakage  in  far  field. 

5.3  Problem  Three  -  Nonlinear  Pressure/Flowrate  in  Far  Field 

Problem  one  and  two  examined  the  extremes  in  far  field  behavior  - 
zero  flow  or  zero  pressure.  In  reality,  it  is  expected  that 
some  leakage  occurs,  but  that  the  leakage  is  pressure  dependant.  We 
examined  this  situation  by  specifying  a  nonlinear  pressure/flowrate 
relation  at  the  far  field  boundaries  (Figure  5.14).  This  relation 
resulted  in  leakage  flow  at  about  one  third  of  the  far  field  nodes. 

Flowrates  at  the  input  node  and  output  nodes  are  shown  in 
Figure   5.15.    The  input  flowrate  is  slightly  large  than  the  flowrate 
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for  problem  one,  which  had  no  leakage.   The  output  flowrate  is  about 

3/4  of  the   input,   so  about  one  fourth  of  the  flow  leaks  to  the  far 

field.  The  pressure  at  a  far  field  node  which  has  no  flow  leakage  is 
shown  in  Figure  5.16. 

Tracer  was  pumped  into  the  crack  system  from  time  1x10   sees  to 

9x10  sees .  The  motion  of  the  tracer  is  shown  in  Figure  5.17. 
Diffusion  of  the  pulse  occurs  between  the  input  and  output.  We 
believe  that  most  of  this  is  realistic  due  to  mixing  of  tracer  in  the 
open  joints.  In  addition,  the  multiple  flow  paths  effectively  diffuse 
the  tracer  pulse. 

Kinetic  energy  during  convergence  is  displayed  in  Figure  5.18. 
Figure  5.19  and  Figure  5.20  show  the  pressures  plot  and  flowrates. 
Figure  5.21  shows  the  velocity  plot.  This  clearly  indicates  the 
slower  flow  velocities  in  the  larger  joints.  In  Fig. 5.20  and 
Fig. 5. 21,  the  arrows  of  fluid  flow  in  far  field  show  1/3  of  far  field 
nodes  leak  due  to  the  nonlinear  relation  of  pressure  vs  flowrate. 

The  tracer  motion  is  shown  in  Figure  5.22  at  three  different  time 
steps.  The  tracer  spread  into  more  fluid  elements  through  time 
history.    Because  we   only  pumped  the  tracer  into  the  well  from  time 

0x10  to  9x10  sec,  Figure  5.22(c)  shows  that  at  time  15x10  sec  the 
zero  concentration  of  tracer  has  spread  from  the  input  area. 

5.4  Problem  Four  -  Reduced  Resistance  Flow  Path 

In  problem  three,  we  examined  flow  for  a  perfectly  symmetric 
situation.  In  reality,  one  might  expect  that  there  will  be  preferred 
flow  paths   of  lower  resistance  between  the  wells,  in  addition  to  the 
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more  uniform  flow  paths.  We  examined  this  situation  in  problem  four 
by  adding  a  third  material  with  an  gap  opening  of  0.0000792  m.  Eight 
elements  which  previously  were  material  type  one  were  changed  to 
material  type  three.  Material  type  three  has  twice  the  gap  opening  of 
material  type  one  (Figure  5.23).  These  eight  elements  construct  a 
"bridge"  between  input  node  and  output  node  which  allows  fluid  flow 
pass  more  easily. 

Figures  5.24(a)  and  (b)  show  the  input  flowrate  and  output 
flowrate.  The  flow  is  increased  over  problem  three  due  to  the 
decreased  resistance  of  the  bridge.  Because  of  the  nonlinear 
pressure  vs  flowrate  relation  in  far  field,  the  leakage  of  fluid  flow 
causes  the  output  flowrate  to  be  less  than  the  input  flowrate. 

The  motion  of  tracer  is  shown  in  Figure  5.25.  The  existence  of 
the  "bridge"  causes  the  tracer  to  move  faster  from  input  node  to 
output  node  through  the  "bridge".  Due  to  the  dilution  of  fluid,  the 
concentration  of  tracer  at  the  output  place  is  very  low  at  the 
beginning  and  then  gradually  increases.  Figure  5.26  shows  the  kinetic 
energy  plot. 

Figure  5.27  shows  the  pressure  plot  and  Figure  5.28  shows  the 
flowrate  plot.  The  velocity  plot  is  shown  in  Figure  5.29.  Figure 
5.28  and  5.29  show  clearly  the  low  resistance  bridge  between  input 
node  and  output  nodes . 

The  motion  of  the  tracer  is  shown  in  Figure  5.30  at  three 
different  time  steps.  The  tracer  follows  the  low  resistance  flow  path 
from  input  node  to  output  node  through  time  history. 
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5.5  Problem  Five  -  Initially  Empty  Fractures 

In  this  problem,  we  used  a  spacing  of  40  m  horizontally  and  40  m 
vertically.  The  larger  joints  are  assumed  to  have  an  opening  of  0.008 
m  (Figure  5.32(a))  and  the  vertical  joints  have  an  opening  of  0.00005 
m  (Figure  5.32(b)).  These  values  were  chosen  to  give  approximately 
the  correct  storage  capacity  and  approximately  the  observed  flowrate. 
It  is  assumed  that  elements  of  material  type  one  (Figure  5.32(a))  are 
active  and  that  of  material  type  two  (Figure  5.32(b))  are  nonactive 
initially,  except  for  those  elements  which  connect  to  the  input  and 
output  nodes  which  are  assumed  to  be  active. 

One  time  step  is  specified  as  40000  seconds  (11.1  hours)  which 
is  less  than  the  time  that  need  to  fill  a  empty  element.  According  to 
the  result,  it  took  640000  seconds  (7.5  days)  for  fluid  to  travel  from 
the  input  node  to  output  node.  This  clearly  indicates  that  the 
assumed  gap  of  the  open  joints  was  too  large. 

Figure  5.33(a)  and  (b)  show  the  input  flowrate  and  output 
flowrate.  The  output  flowrate  is  less  than  the  input  flowrate  because 
of  the  far  field  leakage. 

Figure  5.34  through  Figure  5.36  give  pressure  and  flowrate  plots 
at  three  different  time  steps.  At  each  time  step,  we  see  that  more 
fluid  elements  fill  with  fluid  and  become  active.  Because  of  the 
nonlinear  pressure  vs  flowrate  relation,  some  far  field  nodes  leak,  as 
shown  in  Figure  5.36(b). 
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Figure  5.1:  Mesh  plot  of  application  problem  one 
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Figure  5.2(a):  Material  property  of  material  type  one 
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Figure   5.3(a):    Flowrate   at   input  node 
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Figure   5.3(b):    Flowrate   at   output  node 
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Figure  5.4:  The  pressure  of  a  far  field  node  No. 30 


a   eee  z   eee  4   eee  s    eee  8   eee 

NO     OF     ITERATE       <»„  3, 


Figure   5.5:    Kinetic   Energy  vs   No   of   Iteration  plot 
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Figure  5.6:  Pressure  plot  of  application  problem  one 
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Figure  5.7:  Flowrate  plot  of  application  problem  one 
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Figure  5.8:  Far  field  condition  of  application  problem  two 
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Figure   5.9(a):    The   input   flowrate   of  application  problem  two 
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Figure  5.9(b):  The  output  flowrate  of  application  problem  two 
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Figure   5.10(a):    The  pressure  at  a  far  field  node  No. 30 
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Figure   5.10(b):    The   flowrate   at  a  far   field  node  No. 30 
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Figure  5.11:  Kinetic  Energy  vs  No  of  Iteration  plot 
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Figure  5.12:  Pressure  plot  of  application  problem  two 


Figure  5.13:  Flowrate  plot  of  application  problem  two 
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Figure  5.14:  Far  field  condition  of  application  problem  three 
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Figure  5.15(a):  Flowrate  at  input  node  (problem  three) 


Figure  5.15(b):  Flowrate  at  output  node  (problem  three) 
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Figure  5.16:  Pressure  at  a  far  field  node  No. 30 
which  has  no  flow  leakage 
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Figure  5.17(a):  Tracer  at  input  node  (problem  three) 
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Figure  5.17(b):  Tracer  at  middle  of  crack  (problem  three) 
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Figure  5.17(c):  Tracer  at  output  node  (problem  three) 
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Figure   5.18:    Kinetic   Energy  vs   No   of  Iteration  plot 
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Figure  5.19:  Pressure  plot  of  application  problem  three 


Figure  5.20:  Flowrate  plot  of  application  problem  three 
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Figure  5.21:  Velocity  plot  of  application  problem  three 


Figure  5.22(a):  Tracer  plot  at  4x10   seconds 
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Figure  5.22(b):  Tracer  plot  at  9x10  seconds 


Figure  5.22(c):  Tracer  plot  at  1.5x10   seconds 
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Figure  5.23:  Material  property  of  material  type  three  (problem  four) 
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Figure   5.24(a):    Input   flowrate   of  application  problem  four 
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Figure  5.24(b):  Output  flowrate  of  application  problem  four 
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Figure   5.25(a):    Tracer  at   input  node    (problem  four) 
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Figure   5.25(b):    Tracer   at  middle   of  crack   (problem   four) 
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Figure  5.25(c):  Tracer  at  output  place  (problem  four) 
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Figure   5.26:    Kinetic   energy  vs   No   of   iteration  plot 
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Figure  5.27:  Pressure  plot  of  application  problem  four 


Figure  5.28:  Flowrate  plot  of  application  problem  four 
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Figure  5.29:  Velocity  plot  of  application  problem  four 


Figure  5.30(a):  Tracer  plot  at  1x10   seconds 
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Figure  5.30(b):  Tracer  plot  at  4x10   seconds 
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Figure  5.30(c):  Tracer  plot  at  6x10   seconds 
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Figure  5.31:  Application  problem  five 
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Figure  5.32(a):  Material  property  of  material  type  one 


Figure  5.32(b):  Material  property  of  material  type  two 
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Figure   5.33(a):    Input   flowrate   of  application  problem  five 
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Figure   5.33(b):    Output   flowrate   of  application  problem  five 


91 


Figure  5.34(a):  Pressure  plot  at  4x10   seconds 


Figure  5.34(b):  FIcv/rate  plot  at  4x10  seconds 
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Figure  5.35(a):  Pressure  plot  at  6x10  seconds 


Figure  5.35(b):  Flowrate  plot  at  6x10   seconds 
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Figure  5.36(a):  Pressure  plot  at  1.2x10  seconds 


Figure  5.36(b):  Flowrate  plot  at  1.2x10   seconds 


94 


CHAPTER  VI 


CONCLUSIONS 


In  this  thesis,  we  have  developed  a  finite  element  model  of 
fluid  flow  in  fractured  rock.  The  flow  paths  in  the  rock  are 
modeled  discretely  using  one  dimensional  finite  elements  combined  in  a 
network.  Special  features  include  the  capability  to  model  the 
filling  of  empty  joints,  simulate  far  field  leakage  using  a  nonlinear 
pressure/flowrate  relation,  and  model  the  injection  of  a  tracer  and 
its  distribution  through  the  system.  Interactive  computer  graphics 
allows  the  user  to  easily  specify  the  problem  and  review  the 
results.  Dynamic  relaxation  is  used  to  obtain  a  solution.  This 
robust  scheme  will  be  used  when  the  fluid  model  is  coupled  to  a 
deformable  rock  model  and  the  solution  becomes  nonlinear. 

The  verification  examples  presented  in  chapter  IV  are  idealistic. 
They  demonstrate  most  features  of  the  model. 

The  application  problems  show  the  pressure  and  fluid  flowrate 
distributions  for  fluid  being  pumped  in  one  well  at  high  pressure  and 
removed  from  another  well  at  low  pressure.  Application  problems  one, 
two,  and  three  show  the  significance  of  the  far  field  boundary 
conditions.  Significant  leakage  occurs  if  the  zero  pressure  condition 
is  specified. 

In  problem  four  we  examined  the  effect  of  a  lower  resistance 
bridge  between  the  input  and  output  wells.  This  is  an  attempt  to 
simulate  a  more  realistic  situation,  where  there  is  a  dominant  flow 
path  between  wells   in  addition  to  many  other  flow  paths.   Flow  did 
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indeed  follow  the  lower  resistance  path,  as  shown  by  both  the  flowrate 
and  tracer  calculations.  The  tracer  signal  was  received  at  the  output 
much  sooner.  Smoothing  of  the  signal  occurs  both  because  of  dilution 
with  water  in  the  fractures  and  due  to  mixing  from  the  different  flow 
paths . 

In  problem  five  we  demonstrated  the  filling  of  open  joints.  The 
analysis  showed  the  progressive  filling  and  activation  of  elements 
between  the  input  and  output. 

All  the  application  calculations  were  approximate  solutions  in 
support  of  the  Hot  Dry  Rock  geothermal  energy  project.  The  accuracy 
of  prediction  depends  largely  upon  the  amount  of  correct  geological 
information  available.    The  controlling  factor  of  fluid  flow  is  the 

3 
magnitude   of  the  aperture,   and  since   flowrate  depends  on  (a)  ,  a 

slight  change  in  aperture  can  easily  dominate  any  other  change  in  the 

geometry  of  the   flow  field.    In  general   the  information  such  as 

fracture  aperature,  wall  roughness  and  far  field  conditions  are  little 

known. 

The   finite   element  model  developed  in  this   study   is  quite 

versatile,   and  can  be  coupled  to  the  deformation  of  the  rock  blocks 

easily.   This  coupled  code  will  make  it  possible  to  predict  the  manner 

in  which  joint  openings  change  during  Hot  Dry  Rock  experiments.   After 

coupling   the   fluid  and  structural  models,  the  next  step  will  be  to 

couple  heat   transfer  to  the  solution.   The  final  model  will  aid  the 

engineer   in  developing  a  clearer  picture  of  flow  behavior  in  the  HDR 

reservoir. 
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ABSTRACT 

In  this  thesis,  we  develop  a  finite  element  model  of  flow  in 
fractured  rock  masses.  The  joints  are  modeled  discretely  using  one 
dimensional  finite  elements  connected  in  a  network.  Because  of  the 
robust  solution  scheme,  complicated  nonlinear  joint  networks  can  be 
solved.  The  model  includes  the  filling  of  empty  joints  as  fluid  is 
pumped  through  the  fracture  network.  If  a  joint  is  not  filled,  the 
pressure  in  the  joint  is  assumed  zero  and  the  net  flow  into  the  joint 
is  calculated.  When  a  joint  fills,  the  element  becomes  active  and  is 
included  in  the  flow  calculation.  Far  field  leakage  is  simulated  using 
a  nonlinear  pressure -flowrate  relation.  The  user  can  also  specify 
arbitrary  tracer  input  and  monitor  the  distribution  of  the  tracer 
during  an  analysis. 

Example  calculations  include  solutions  in  support  of  the  Hot  Dry 
Rock  geo thermal  energy  project.  The  solutions  show  the  pressure 
and  fluid  velocity  distributions  for  fluid  being  pumped  in  one  well 
at  high  pressure  and  removed  from  another  well  at  low  pressure.  A 
comparison  of  solutions  with  uniform  and  nonuniform  joint  openings 
shows  how  the  fluid  follows  the  low  resistance  path  and  how  the 
tracer  output  signal  is  smoothed  as  a  result  of  the  different  flow 
path  arrival  times.  The  results  of  this  model  will  aid  engineers  of 
Los  Alamos  to  understand  more  about  fluid  motion  in  the  Hot  Dry  Rock 
reservoir. 


